Transmission and Drug Resistance Genotype of Multidrug-Resistant or Rifampicin-Resistant Mycobacterium tuberculosis in Chongqing, China

ABSTRACT Multidrug-resistant or rifampicin-resistant tuberculosis (MDR/RR-TB) is a global barrier for the Stop TB plan. To identify risk factors for treatment outcome and cluster transmission of MDR/RR-TB, whole-genome sequencing (WGS) data of isolates from patients of the Chongqing Tuberculosis Control Institute were used for phylogenetic classifications, resistance predictions, and cluster analysis. A total of 223 MDR/RR-TB cases were recorded between 1 January 2018 and 31 December 2020. Elderly patients and those with lung cavitation are at increased risk of death due to MDR/RR-TB. A total of 187 MDR/RR strains were obtained from WGS data; 152 were classified as lineage 2 strains. Eighty (42.8%) strains differing by a distance of 12 or fewer single nucleotide polymorphisms were classified as 20 genomic clusters, indicating recent transmission. Patients infected with lineage 2 strains or those with occupations listed as “other” are significantly associated with a transmission cluster of MDR/RR-TB. Analysis of resistant mutations against first-line tuberculosis drugs found that 76 (95.0%) of all 80 strains had the same mutations within each cluster. A total of 55.0% (44 of 80) of the MDR/RR-TB strains accumulated additional drug resistance mutations along the transmission chain, especially against fluoroquinolones (63.6% [28 of 44]). Recent transmission of MDR/RR strains is driving the MDR/RR-TB epidemics, leading to the accumulation of more serious resistance along the transmission chains. IMPORTANCE The drug resistance molecular characteristics of MDR/RR-TB were elucidated by genome-wide analysis, and risk factors for death by MDR/RR-TB were identified in combination with patient information. Cluster characteristics of MDR/RR-TB in the region were analyzed by genome-wide analysis, and risk factors for cluster transmission (recent transmission) were analyzed. These analyses provide reference for the prevention and treatment of MDR/RR-TB in Chongqing.

tried to improve surveillance of recent transmission and drug genotypic resistance patterns of MDR/RR-TB using WGS in Chongqing. Our research may promote MDR/RR-TB control in Chongqing and other similar regions.

RESULTS
Treatment results of MDR/RR-TB patients. From 1 January 2018 to 31 December 2020, we diagnosed 223 (24.6%) MDR/RR-TB patients by phenotypic drug susceptibility testing (DST) from a total 906 culture-positive TB patients in the Chongqing Tuberculosis Control Institute. Treatment information for 200 MDR/RR-TB patients with culture-positive isolates was collected. As shown in Table 1, except for 2 patients who died from other causes, 23 (11.5%) patients died from MDR/RR-TB, and 80 (40%) patients received favorable outcomes (34 cured and 46 completed the treatment). Ninety-five (47.5%) patients were classified into the pending group for various reasons, with 49 (24.5%) still undergoing treatment, while the reasons for the other 46 cases include loss to follow-up, severe adverse drug reactions, and treatment refusal.
Sociodemographic characteristics (gender, age, residence, occupation, etc.) and clinical characteristics (complications, previous TB treatment history etc.) of patients in those groups from Table 1 (excluding the 2 deaths from other causes) are shown in Table 2. Univariate and logistic regression analyses were also performed for these 198 patients. As shown in Table 3, the results indicated that age, occupation, and lung cavitation were significant between favorable outcome and death (P # 0.1).
Next, these variables were progressed to the multivariate analysis (Table 4). Age was the risk factor for death for elderly patients. Between the pending and death groups, the results from the univariate analysis indicated that age, occupation, and residence were significant risk factors (P # 0.1). These variables were progressed to the multivariate analysis ( Table 4). The age was the risk factor for older patients. Between the favorable outcome and pending groups, the results from the univariate analysis indicated that age, residence, and history of treatments were significant (P # 0.1). Then, these variables were progressed to the multivariate analysis ( Table 4). History of treatments was the risk factor; the patients previously treated with second-line drugs were more likely to have pending outcome.
Recent transmissions of MDR/RR-TB in Chongqing, China. Incident TB is a combination of recent transmission and remote transmission (initial transmissions of a non-MDR/RR strain that developed into MDR/RR-TB over a period of years) (21). To evaluate the burden of recent transmission of MDR/RR-TB in Chongqing, a minimum distance spanning tree (MST) based on SNP distance was calculated on the basis of the WGS sequencing data. A total of 187 isolates were obtained from 200 culture-positive MDR/ RR isolates. Eighty (42.8%) isolates differed in 20 genomic clusters by an allele cutoff of 12 or fewer based on the cg-MLST approach (see Fig. S1 in the supplemental material). The clusters, ranging from 2 patients (12 clusters) to 16 patients (1 cluster), presented with recent transmission (Fig. 1). Spatial-location analysis revealed that 43 (53.8%) clustered isolates were located in the same district (Fig. 2). Pearson correlation analysis demonstrated minimal association between spatial distance and genetic distance (SNPs) for all 80 genomic clustered strains (CRR = 20.09; P = 0.14). Only in cluster 4, a strong positive correlation (CRR) was observed between spatial distance and genetic distance (SNPs) (CRR = 0.95; P , 0.01). (Details are shown in File S1.) Of all the MDR/RR isolates, 152 (85.4%) were grouped into lineage 2. As shown in Fig. 1, within the 80 isolates in 20 genomic clusters, 74 isolates (17 clusters) were lineage 2 and 6 isolates (3 clusters) were lineage 4. Then, we analyzed the risk factors for recent transmission. The results of univariate analysis suggested that gender, occupation, residence, and lineage were significantly associated with disease transmission (P # 0.1) (Table 5). Then, gender, occupation, residence, and lineage were progressed to multivariate analysis (Table 6). Occupation and lineage were the risk factors by which the patients with other occupations (mainly household management service industry) or infected by lineage 2 strains were more likely to be part of a cluster. The risk factors of patients who were infected by lineage 2 were analyzed (data not shown).  Cluster transmission was the only risk factor for Mycobacterium tuberculosis lineage 2 infection. In conclusion, cluster transmission of lineage 2 strains by patients with other occupations is a significant risk factor driving the epidemic toward MDR/RR-TB. Drug resistance profile of MDR/RR-TB. For all 200 culture-positive MDR/RR isolates, we calculated the rates of resistance to 12 antibiotics ( Fig. 3; MIC details are in File S3). One hundred (50.0%) strains were MDR strains that were resistant to isoniazid (INH) and rifampicin (RIF). Ninety-one (45.5%) MDR/RR-TB cases were classified as pre-extensively drug resistant (pre-XDR), which are additionally resistant to one of the fluoroquinolones (FQs) (22). One hundred thirty-two (66.0%) strains were resistant to ethambutol (EMB), which is one of four first-line antituberculosis drugs. The rates of resistance to two FQs, LVX and MXF, were 45.5% and 48.0%. The rates of resistance to two aminoglycosides (AGs), kanamycin (KAN) and amikacin (AMK), were 14.5% and 13.5%. For new antituberculosis drugs, the rates of resistance to BDQ, delamanid (DLM), LZD, and CFZ were 1.5%, 4.5%, 1.5%, and 0.5%, respectively. Five (2.5%) MDR/RR-TB cases were classified as extensively drug resistant tuberculosis (XDR-TB), which, by definition, is resistant to a fluoroquinolone and either BDQ or LZD (or both) (22). Lineage 2 was more likely to be resistant to ethambutol than lineage 4 ( Table 7).
Transmission networks of MDR/RR-TB. A phylogenetic tree of 187 MDR/RR-TB isolates indicating drug resistance profiles and lineages is shown in Fig. 4. A total of 152 isolates belonged to lineage 2 (East Asian), and the rest belonged to lineage 4 (Euro-American). As illustrated in Table 5, the isolates from lineage 2 were the main strains that resulted in the recent transmission of MDR/RR-TB in Chongqing.
We observed the interesting phenomenon that isolates with the same rpoB mutations were scattered in clusters in the phylogenetic tree (Fig. 4). To graphically certify the transmission of an MDR/RR strain, we constructed a median-joining network of each genomic cluster and mapped genotypic drug resistance onto the network (Fig. 5). The strains in 16 genomic clusters (including 2 clusters in which the strains had no rpoB mutations) had the same or no RIF resistance-associated rpoB mutations individually, except cluster 3 (one strain had rpoB_p.His445Ser plus rpoB_p.Leu430Pro, while others obtained rpoB_p.Ser450Leu), cluster 5 (one strain contained rpoB_p.Asp435Val, while others obtained rpoB_p.Ser450Leu), cluster 14 (the two strains had rpoB_p.Asp435Val and rpoB_p.His445Leu individually), and cluster 6 (each of 2 strains had the same mutation sites). Therefore, 76 (95.0%) of all 80 strains defined as recent transmission have the same or no mutations in the rpoB gene within each cluster. For INH, at least two strains shared   (Fig. 5).

DISCUSSION
In Chongqing, 23 of all 200 MDR/RR-TB patients died. The risk factor analysis for death showed that age was the major risk factor. Consistent with our research, MDR-TB patients aged 60 and over or with cavitary disease are more likely to have poor treatment outcomes in Zhejiang, China (23). Treatment history was the risk factor for the patients previously treated with firstor second-line drugs, and these patients were more likely to have pending outcome. Therefore, elderly patients or those previously treated with the firstor second-line drugs should have enhanced management and supervision to complete standard treatment regimens. The control of TB transmission depends on the identification and treatment of infectious patients and their close contacts (24). In China, a setting with a high burden of TB, epidemiological data are often difficult to get. WGS can identify reasonable transmission events in patients without prior recourse to epidemiological data (5,14,25,26). Among WGS analysis pipelines for detection of epidemiologically linked tuberculosis cases, the cg-MLST approach is standardized and easy to get compared to the SNP-based pipelines (27). In Chongqing, 80 (42.8%) of all 187 MDR/RR-TB patients were identified by recent    (13). Another study reported that the local incidence of TB in urban centers occurs via local transmission between both migrants and residents in Shanghai (28). In India, transmission of particular pre-XDR/XDR lineage 2 strains is the main driving force of the pre-XDR/XDR-TB epidemic (29). In our research, the diagnosis delay defined by the interval between the diagnosis of TB and the diagno- Furthermore, another study demonstrated that WGS can predict susceptibility of Mycobacterium tuberculosis to first-line drugs more accurately than phenotypic testing (15). For the new and repurposed tuberculosis drugs, no resistance mutations were observed in phenotypically BDQ-, DLM-, LZD-, and CFZ-resistant strains. The mutations associated with resistance to the new and repurposed Mycobacterium tuberculosis drugs should be updated in their resistance gene catalogue database (31)(32)(33)(34).
The genomic clusters might represent transmission of an MDR/RR strain or initial transmission of a non-MDR/RR strain that later developed multidrug resistance (13). It is known that drug resistance of tuberculosis is mainly conferred by mutations in genes coding for drug targets or converting enzymes (35). Prior mutations may accumulate along the transmission chain, and strains may acquire new mutations that further increase drug resistance. An earlier report highlighted that both drug resistance transmission and amplification contribute to disease burden globally (36). In China, it has been proved that additional drug resistance amplified for Mycobacterium tuberculosis during the turnaround time for drug susceptibility testing (37). In Chongqing, 55% of the MDR/RR-TB strains accumulated additional drug resistance mutations along the transmission chain, especially against FQs (63.6% [28 of 44]). It is suggested that most strains developed additional drug resistance during the expansion of MDR/RR clusters, indicating that MDR/RR TB is not fully controlled. The use of a WGS-based approach for surveillance purposes might enable the public health service to take appropriate control actions in specific settings in China.

MATERIALS AND METHODS
Strain collection. All suspected pulmonary tuberculosis patients are referred to local appointed hospitals, where the diagnosis is made by sputum smear and culture. For isolation from culture, each specimen was treated with 1 volume of 4% sodium hydroxide per 1 volume of sputum and then homogenized by vigorous stirring. An aliquot of 0.1 mL of the resulting specimen was inoculated into two tubes of acidified Löwenstein-Jensen (L-J) medium and incubated at 37°C. The culture was first assessed during week 1 for rapidly growing bacteria and every week for slowly growing bacteria; if no bacterial growth was observed by week 8, the result was recorded as negative. All rifampicin-resistant Mycobacterium tuberculosis isolates were previously identified using the proportion method on L-J medium containing rifampicin at a concentration of 40 mg/mL. All cases of MDR/RR-TB are reported to the Chongqing Tuberculosis Control Institute, and all culture-positive samples are delivered to the national reference laboratory for tuberculosis (NTRL) in the China CDC.
From 1 January 2018 to 31 December 2020, we enrolled 223 (24.6%) MDR/RR-TB patients from a total of 906 culture-positive TB patients in the Chongqing Tuberculosis Control Institute. The strains were thawed and subcultured on L-J medium for further analysis by combining phenotypic drug susceptibility testing (DST) and WGS at the NTRL. A total of 200 MDR/RR strains isolated from unique TB patients were used in our population-based, retrospective observational study. The strains with reculture failure and serial samples from identical patients were excluded from this study.
Patient information collection. Sociodemographic characteristics (gender, age, residence, and occupation) and clinical characteristics (complications and previous TB treatment history) of patients were extracted from the national drug resistance surveillance database; this information was collected and compiled by local physicians using a questionnaire filled in by patients after written informed consent was obtained at the time of patient visits.
Strain identification. Matrix-assisted laser desorption ionization-time of flight mass spectrometry (MADLI-TOF MS) was used to distinguish Mycobacterium tuberculosis complex (MTBC) from other mycobacteria. The detailed protocol is as follows.
Ten microliters of standard ring bacteria was harvested from L-J medium and dispersed in 1 mL of 75% alcohol (high-performance liquid chromatography [HPLC] grade), mixed well, and stored at 220°C for use. Prior to analysis, the bacterial suspension was centrifuged and the sediment was resuspended in 10 mL of acetonitrile with a small amount of zirconia/silica microbeads. Next, 10 mL of 70% formic acid was added after full-speed vortexing. The supernatant was reserved after centrifugation.
For each sample, 1 mL was deposited on a polished steel MSP 384 target plate (Bruker Daltonics, Bremen, Germany) and 1 mL of matrix solution (saturated a-cyano-4-hydroxycinnamic acid [HCCA]) was then added. The samples were air dried for 5 min before being processed in the mass spectrometer. To validate the analysis of a whole MSP 384 target, bacterial test standard (Escherichia coli DH5a protein extract) was used as a positive control and noninoculated matrix solution (HCCA) was used as a negative control. The analyses were performed using flexControl 3.0 software (Bruker Daltonics). The spectra were analyzed within an m/z range of 2,000 to 20,000. Four raw spectra were automatically acquired using the flexControl 3.0 software and then compared with the Bruker Daltonics database using MALDI Biotyper 3.0 software. To validate the analysis using MALDI Biotyper software, the identification of the positive control was required to be E. coli with an identification score of $2, and the negative control had to yield a nonidentifying score of #1.7.
Risk factors of death. In our research, all patients were classified into four groups based on the treatment outcomes when we collected the information. In the favored outcome group, the patients were cured (defined as a bacteriologically confirmed TB patient who was smear or culture negative in the last month of treatment and on at least one previous occasion during treatment) or completed treatment (defined as a bacteriologically confirmed TB patient who completed treatment without evidence of failure but with no record of a negative sputum smear or culture from the last month of treatment and at least one previous occasion during treatment); in the pending group, the patients had no treatment outcome (patients were in treatment, had an adverse reaction, refused treatment, or were lost to follow-up; the other two groups consisted of patients who died from TB and those who died from other causes. Univariate and logistic regression analyses were performed for 198 patients between each group, except for the group of death by other causes (2 patients). When presenting a risk factor analysis, we first showed the results from the univariate analysis, indicating which variables were significant (P , 0.1). Then, these variables were progressed to the multivariate analysis, and then we showed the results for the multivariate analysis (P value of ,0.05 was significant).
DST. Drug susceptibility testing (DST) of Mycobacterium tuberculosis strains against rifampicin, isoniazid, ethambutol, kanamycin, amikacin, moxifloxacin, levofloxacin, ethionamide, bedaquiline, delamanid, linezolid, and clofazimine was performed using a UKMYC5 plate (Thermo Fisher Scientific Inc., USA), which has been reported as an alternative DST method with high accuracy and reproducibility (38). DST was conducted strictly according to the manufacturer's instructions by trained staff at the national tuberculosis reference laboratory of China.
Definitions. MDR/RR-TB was defined as Mycobacterium tuberculosis resistance to at least rifampicin. MDR-TB was defined as Mycobacterium tuberculosis resistance to at least isoniazid and rifampicin. MDR-TB only was defined as an MDR-TB strain that was susceptible to both fluoroquinolones (moxifloxacin or levofloxacin) and the second-line anti-TB drugs (amikacin or kanamycin). Pre-XDR-TB was defined as MDR-TB with additional resistance to any fluoroquinolones (moxifloxacin or levofloxacin). XDR-TB was defined as MDR-TB with additional resistance to any fluoroquinolones and either BDQ or LZD.
DNA extraction and sequencing. We used the cetyltrimethylammonium bromide (CTAB) method to extract genomic DNA from cultures of one sputum specimen per patient. Sequencing libraries were prepared by using the Illumina Nextera kit following the manufacturer's protocol and sequenced on the Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) with 2 Â 150 paired-end (PE) strategies. Coverage of 100Â was expected. All whole-genome sequencing procedures were performed by Annoroad Gene Technology Company (Beijing, China). In total, we obtained WGS data for 187 strains.
Bioinformatics analysis. To guarantee good quality of the sequencing read, specific parameters were set and followed: 60 to 65% GC content of original data, #20% per-base sequence content, #2% overrepresented sequences, #10% reads containing joint sequences, #2% low quality reads, $95% of the reads mapping to reference genome, #20% duplicate reads, $95% of the reference genome being mapped by reads, $50Â average genome sequencing depth, and #1% base mismatch for the reads mapping to reference genome.
We used the functionality implemented in the Ridom SeqSphere1 software (version 7.2.3) with default settings to perform cg-MLST analysis. The genome of the Mycobacterium tuberculosis strain H37Rv (GenBank accession no. NC_00962.3) served as the reference genome. Afterwards, gained genomes were compared to the seed genome to identify a list of core genome genes. Here, default settings include the removal of the shorter of two genes overlapping by more than four bases and of genes with an internal stop codon in .90% of all genomes from the scheme. cg-MLST-based neighbor joining and minimum spanning trees (MST) were calculated and drawn with SeqSphere1 software.
Risk factors of cluster transmission (recent transmission). Cluster transmissions were defined as the isolates with pairwise genetic distances of fewer than 12 SNPs. Univariate and logistic regression analyses were performed for 187 patients between cluster group and noncluster group. During the presentation of a risk factor analysis, we first showed the results from the univariate analysis, indicating which variables were significant (P , 0.1). Then, these variables were progressed to the multivariate analysis, and then we showed the results for the multivariate analysis (P value of ,0.05 was significant).
Spatial-location analysis. The residential address listed in the questionnaire of each patient was geocoded with QGIS 3.12 and Baidu Maps (Baidu, Beijing, China) to verify locations. Spatial distance indicated the distance between where two patients live on the map. Genetic distance indicated the number of different SNPs between two strains identified by the cg-MLST approach. Pearson correlation analysis between spatial distance and genetic distance (SNPs) for all 80 genomic clustered strains was performed by IBM SPSS statistics software. The greater the absolute value (0 to 1) of correlation (CRR), the stronger the correlation. Positive value is positive correlation, and negative value is negative correlation. A P value of ,0.01 indicated significant correlation.
Statistical analysis. Descriptive statistics was performed for patients' demographics and lineages, resistance categories, and clustering status of MTBC strains. Data stemming from genomic analysis of clinical isolates were analyzed statistically using IBM SPSS statistics software for Windows (version 19) and R (version 3.6.1). For univariate analysis of potential factors, we performed Fisher's exact test. Factors of a significant result from the univariate model were included into a multivariate logistic regression analysis. Odd ratios with 95% confidence interval (CI) were estimated, and variables with P values of less than 0.05 were taken as significant predictors.
Ethics statement. The institutional review boards of the China CDC approved the study. All patients provided written informed consent.
Data availability. The accession numbers of all sequenced genomes from this study are shown in File S4. All the staff from National Tuberculosis Reference Laboratory in China is also highly acknowledged. We thank Ruichao Yue and Wei Guo from the University of North Carolina at Greensboro for the revision of the language.